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Abstract 

In mathematical finance a popular approach for pricing options under some Levy model is to 
consider underlying that follows a Poisson jump diffusion process. As it is well known this results 
in a partial integro-differential equation (PIDE) that usually does not allow an analytical solution 
while numerical solution brings some problems. In this paper we elaborate a new approach on how 
to transform the PIDE to some class of so-called pseudo-parabolic equations which are known in 
mathematics but are relatively new for mathematical finance. As an example we discuss several 
jump-diffusion models which Levy measure allows such a transformation. 
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1 Introduction 



In mathematical finance a popular approach for pricing options under some Levy model is to consider 
underlying that follows a Poisson jump diffusion process. As it is well known this results in a partial 
integro-differential equation (PIDE) that usually does not allow an analytical solution while numerical 
solution brings some problems. These problems are mainly related to computing a non-local integral 
term while we assume that computing a differential part of the PIDE, being discussed numerous 
times in the literature, could be provided in a relatively standard way. Moreover, using splitting 
technique it is always possible to reduce the whole PIDE to a series of equations part of which are 
pure PDE and the rem a ining p art are pure evo l utiona ry-integral equations (EIDE) (see, for instance, 
in 't Hout and Welfert ( 20091 ). It kin and Carr ( 20061 )). Thus, further on we will consider jus t the 



(2009), 


Hilber et al. 


(2009h while i 


in 


Carr and Mavo (20071). Strauss! 



According to the last cited paper we could distinguish the following methods that were used to 
solve the EIDE. In an early paper, Amim (1993) used an explicit multinomial tree based approach. 
DHalluin et al. (2004, 2005b) implemented implicit methods for evaluating vanilla European options, 
barrier options, and American options. They also showed that when a log spaced grid is used with 
a Crank Nicolson discretization on a problem with constant parameters the resulting scheme is un- 
conditionally strictly stable. In addition, they showed that the simple Picard iteration scheme (also 
suggested by Tavella and Randall (2000)) for solving the discretized equations is globally convergent. 
Specifically, they reported that when they priced options in the Merton model the error was reduced 
by two orders of magnitude at each iteration for typical values of the time step size and Poisson arrival 
intensity. More recently, dHalluin et al. (2005a) presented a semi-Lagrangian approach for pricing 
American Asian options under jump diffusion processes. Andersen and Andreasen (2000) derived a 
forward equation describing the evolution of European call options as functions of strike and maturity, 
and discussed its application to the problem of fitting the stock process to option prices in the market. 
They also presented a second order accurate unconditionally stable operator splitting (ADI) method 
for pricing options which does not require iterative solution of an algebraic equation at each time step. 
(Unfortunately, it is not clear how to extend their method to the valuation of American options while 
retaining second order accuracy.) Cont and Voltchkova (2003) used a discretization that is implicit 
in the differential terms and implicit in the integral term, and showed that it converges to a viscosity 
solution. Their method extends to infinite activity models, and does not require the diffusion part of 
the equation to be non-degenerate. These partial integro-differential equations have also been solved 
by many others. See, for example, Zhang (1993) and Matache et al. (2002). Although the pricing 
equations have often been solved numerically, because of the integrals in the equations the methods 
have proven relatively expensive. The obvious discretizations of the pricing equations combine stan- 
dard discretization methods for the differential terms with quadrature methods such as Simpsons rule 
or Gaussian quadrature for evaluating the integral term. This approach is computationally expensive 
since the integral must be approximated at each point of the mesh used for discretizing the differential 
terms. The difficulties are greater if an implicit discretization of both the integral and the differential 
terms is used. The expense of evaluating the integral at all points of the computational grid can, 
however, be reduced by making the same exponential change of variables often used when solving the 
BlackScholes differential equation when there is no jump process. This converts the integral term into 
a correlation integral which can be evaluated at all the mesh points simultaneously using the Fast 
Fourier Transform. This approach has been suggested by many authors ( Wilmott, 1998; Tavella and 
Randall, 2000; Andreasen and Anderson, 2000). 

For multidimensional Levy process various kind of finite elements methods were proposed (see 
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survey in iHilber et al.1 (|200gh ) because finite difference methods are not efficient when dimensionality 
of the problem exceeds 3. 

Also note a new method for exponential jumps proposed by Lipton and Sepp ( 20091 ) who calculate 
the jump integral recursively on the spatial grid. This is a special trick for exponential jumps, it 
does not work for more familiar Gaussian jumps. The authors claim that for discrete jumps a simple 

interpolation routine is sufficient. 

As Carr and Mayo mentioned (see Carr and Mayol (j2007l )) quadrature methods are expensive since 
the integrals must be evaluated at every point of the mesh. Though less so, Fourier methods are also 
computationally intensive since in order to avoid wrap around effects they require enlargement of the 
computational domain. They are also slow to converge when the parameters of the jump process 
are not smooth, and for efficiency require uniform meshes. Therefore, they proposed a different and 
more efficient class of methods which are based on the fact that the integrals often satisfy differential 
equations. Depending on the process the asset follows, the equations are either ordinary differential 
equations or parabolic partial differential equations. Both types of equations can be accurately solved 
very rapidly. They used to demonstrate the advantage of such an approach for the Merton and Kou 
models. However, for other types of the Levy models an extension of their idea is unknown yet. 

Therefore in this paper we propose two different approaches. The idea of the first one is to 
represent a Levy measure as the Green's function of some yet unknown differential operator A. If we 
manage to find an explicit form of such an operator then the original PIDE reduces to a new type of 
equation - so-called pseud o -para bolic equation. These equations are known in mathematics (see, for 



instance, Cannon and Lin ( 19881 )) but are new for mathematical finance. 



Then we rely on two important results, namely: a) the inverse operator A' 1 exists, and b) the 
obtained pseudo parabolic equation could be formally solved analytically via a matrix exponent. 
Having that we discuss a numerical method of how to compute this matrix exponent. We show that 
we can do it using a finite difference scheme similar to that used for solving parabolic PDEs and the 
matrix of this FD scheme is banded. We fulfill this program for general tempered stable processes 
(GTSP) with an integer damping exponent a. 

Alternatively for some class of Levy processes, known as GTSP/KoBoL/SSM models, with the real 
dumping exponent a we show how to transform the corresponding PIDE to a f ractional PDE (method 



2). Fr actional PDE s for the Levy processes with finit e varia tion were derived by lBoyarchenko and Levendorskii 



(|2002h and later bv lCartea and del Castillo-Negrete! (jioO^L 

using a characteristic function techn ique- 
Numerical solution of t h ese eq uations was investigated by Cartea and del Castillo-Negrete ( 20071 ) and 
Marom and Momoniat ( 20091 ). In this paper we derive them in all cases including processes with 



infinite variation using a different technique - shift operators. Then to solve them we apply a new 
method, namely: having results computed for a £ I we then interpolate them with the second order 
in a to obtain the solution at any a£l. 

We also show that despite it is a common practice to integrate out all Levy compensators in the 
integral term when one considers jumps with finite activity and finite variation, this breaks the stability 
of the scheme, at least for the fractional PDE. Therefore, in order to construct the unconditionally 
stable scheme one must keep the other terms under the integrals. To resolve this in Cartea (2007) 
the authors were compelled to change their definition of the fractional derivative. 

We also propose the idea of solving FPDE with real a by using interpolation between option prices 
computed for the closest integer values of a. For the latter an efficient scheme is proposed that results 
in LU factorization of the band matrix. 

It is important to note that both proposed methods could be easily generalized for a time- 
dependent Levy density. 

The rest of the paper is organized as follows. In section [2] we discuss a basic example of the 
method which is built based on a simple exponential Levy measure. In the next section we consider 
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GTSP models and show how to reduce the corresponding PIDE to a pseudo parabolic equation in this 
case. Section U] describes numerical solution of the obtained pseudo parabolic equations in case a E I. 
Section [5] describes a general case of real a and introduces our method of deriving fractional PDE based 
on shift operators. In section [6] we discuss how to solve these FPDE by constructing unconditionally 
stable finite difference schemes of high order of accuracy in space and time and provide some numerical 
examples and comparison with the other methods. The last section concludes. 

2 Basic model 

In this section we consider the simplest possible problem to demonstrate basics of our new method. 
We assume no arbitrage so that there exists a risk-neutral measure Q. We assume zero interest rates 
and dividends so that the stock price is a Q-martingale. Suppose that the underlying stock price 
process is pure jump (i.e. there is no continuous martingale component). Further suppose that the 
jump process is a compound Poisson process. The arrival rate of a jump is constant at A > 0, while 
the jump size distribution is a symmetric Laplace distribution, i.e the probability density for a jump 
of size j £ R, given that a jump has occurred is given by: 

q(j) = — j— , j e R, (i) 

where a > is a free parameter. We recognize that these dynamics let prices become negative and 
ignore this complication. Let u(x,t) be the value of the contingent claim at calendar time t £ [0, T] 
given that the time t stock price is x £ R. As a result of our assumptions, the contingent claim value 
solves the following PIDE: 

d f d 

—u(x, t) + A / \u{x + j, t) - u(x, t) - —u(x, t)j]q(j)dj = 0, (2) 
ot J R dx 

on the domain x G R, t £ [0, T]. For a European call, the terminal condition is 

u(x,T) = (x-K) + , xGR, (3) 

where K £ R is the strike price. 

Now the symmetry of the PDF in ([T]) implies that: 



jq(j)dj = 0, (4) 



and hence the PIDE ([2]) simplifies to: 





(x, t) - Xu(x, t)+\ u(x + j, t)q(j)dj = 0. (5) 
it 



If we do the change of variable z = —j in the integral, we obtain a convolution: 





u(x,t) — Xu(x,t) + A / u(x — z, t)q(z)dz = 0. (6) 

at 



If we do the change of variable y = x — z in the integral, we obtain: 





u(x,t) - Xu(x,t) + X j u(y,t)q(x - y)dy = 0. (7) 
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Now consider the simple second order linear inhomogeneous ODE: 

g"(x) - a 2 g(x) = -S(x), x G R, (8) 

where 5(x) denotes Dirac's delta function. Suppose that the ODE is to be solved subject to the 
boundary conditions: 

lim g(x) = 0. (9) 

x— >±oo 

The solution to this problem is usually referred to as a Green's function. The solution is well known 
to be: 

„— a\x\ 

9{x) = (10) 

Comparing (fTOj) and (P), we see that: 

q(x) = a 2 g(x). (11) 

Hence, the PIDE (|16p can be re-written as: 

d f 

—u(x,t)-\u(x,t) + \a 2 u(y,t)g(x - y)dy = 0. (12) 
at Jr 

To exploit the connection (jlip . let T> x denote the first derivative operator and let A x denote the 
following linear differential operator: 

A x = Vl - a 2 X x , (13) 

where I x is the identity operator. 

Using this operator notation, the ODE ([8]) reads: 

A x g(x) = -5(x). (14) 

Suppose that we apply the A x operator to (fT2|) : 

d f 

A x —u(x, t) - XA x u(x, t) + Xa 2 / u(y, t)A x g(x - y)dy = 0. (15) 

where we have assumed that the interchange of the integral and the differential operator is permissible. 
Substituting {H]) in ([12]) implies that: 

d f 

A x —u(x,t)-\A x u(x,t)-Xa 2 u(y,t)8(x - y)dy = 0. (16) 

Using the sifting property of the delta function implies that our problem reduces to a (third order) 
PDE: 

d 

A x —u(x,t) — \A x u(x,t) — Xa 2 u(x,t) = 0. (17) 
Substituting (I13p in ()1T|) and simplifying implies: 

g^2Q- t u(x, t) - a 2 -Q t u(x, t) - A^n(x, t) = 0. (18) 

Note that the generalization from exponential type kernels to Erlang type kernels can be handled 
by replacing the second order differential operator A x by a higher order differential operator. We 
further note that the Central Limit Theorem implies that the limiting sum of these independent 
exponential random variables is normally distributed. A Gaussian component to the jump kernel 
induces an infinite order ODE which is equivalent to a PDE. Hence the Gaussian type jump of 
Merton can be handled by solving a PDE as we already know. The PDF of a linear combination of 
independent exponential and Gaussian random variables is called the Polya Laguerre distribution. A 
good reference for the above inversion is Hirschman and Widder. 
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3 GTSP /KoBoL /SSM model 



Stochastic skew model (SSM) has been proposed bv lCarr and Wu" ( 2004 ) for pricing currency options. 
It makes use of a Levy model also known as generalized tempered stable processes (GTSP) (see 
Cont and Tankov ( 20041 ) ) for the dynamics of stock prices which generalize the CGMY pr o cesses 
proposed by Carr et al.l (2002). A simil a r mod el was independently proposed by Koponen ( 19951 ) 
and then Boyarchenko and Levendorskii ( 20021 ) The processes are obtained by specifying a more 
generalized Levy measure with two additional parameters. These two parameters provide control on 
asymmetry of sm all jumps and different frequencies for upward and downward jumps. The results of 
Zhou et al.l (|2005l ) show that this generalization allows for more accurate pricing of options. 

Generalized Tempered Stable Processes (GTSP) have probability densities symmetric in a neigh- 
borhood of the origin and exponentially decaying in the far tails. After this exponential softening, 
the small jumps keep their initial stable-like behavior, whereas the large jumps become exponentially 
tempered. The Levy measure of GTSP reads 



-v-\y\ 



l+E^J-jKO + A +i |l+a+ X ^>0' 



-v+\y\ 



\y\ \y\ 

where u± > 0, A± > and a± < 2. The last condition is necessary to provide 



(19) 



y 2 fi(dy) < oo, / n(dy) 

■1 J\y\>l 



< oo. 



(20) 



The case A+ = A_,a + = a_ corresponds to the CGMY process. The limiting case a+ 
0, A + = A_ is the special c ase of the Varian ce Gamma process of Madan and Seneta ( 1990l ). As 
Hagan at al mentioned (see IZhou et al.l (|2005l )) six parameters of the model play an important role 
in capturing various aspects of the stochastic process. The parameters A± determine the overall and 
relative frequencies of upward and downward jumps. If we are interested only in jumps larger than 
a given value, these two parameters tell us how often we should expect such events. u± control the 
tail behavior of the Levy measure, and they tell us how far the process may jump. They also lead 
to skewed distributions when they are unequal. In the special case when they are equal, the Levy 
measure is symmetric. Finally, a± are particularly useful for the local behavior of the process. They 
determine whether the process has finite or infinite a ctivity, or variation. 

Using this model of jumps I Carr and Wu. (120041 ) derived the following PIDE which governs an 
arbitrage-free value of a European call option at time t 



r d C(S, V R , V L ,t) = -C(S, V R , V L ,t) + (r d - r f )S—C(S, V R , V L ,t) 



(21) 











+ k(1 - V R )—C(S, V r , V L ,t) + k(1 - V L )—C(S, V R , V L ,t) 



dV R 

^s\Vr + v l ) d 2 
2 ds 2 



dV L 

C(S,VR,V L ,t) + ap R a v SV R 



u 2 



dSdV R 



C(S,V R ,V L ,t) 



+ ap L a v SV L 
+ \/Vr~ 



C(S, Vr, V L ,t) + ^^C(S, V R , V L ,t) + ^L^CiS, Vr, V L ,t) 



2 dVi 



dSdV L 

c(s e y,VR,v L ,t)-c(s,VR,v L ,t) 
c(s e y,VR,v L ,t)-c(s,VR,v L ,t) 



d_ 

dS 
d_ 
dS 



2 dVt 

c(s,v R ,v L ,t)s( e y - 1) 
c(s,v R ,v L ,t)s( e y -i) 



-vr\v\ 



\l+a 



dy 



e -VL\y\ 
A , M , — dy, 
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on the domain S > 0, V R > 0, Vl > and t € [0,T], where <S, Vr, Vj, are state variables (spot price 
and stochastic variances). For the following we make some critical assumptions. 

1. This PIDE could be generalized with allowance for GTSP processes, which means we substitute 
a in Eq. (j2Tj) with a R , ocl-, and A with Ar, Xl correspondingly. 



2. The obtained PIPE c ould be solved by using a splitting technique similar to that proposed in 
Itkin and CaTrl \20O& . 



3. We assume a R < 0, oll < which means we consider only jumps with finite activity. Therefore, 
each compensator under the integral could be integrated out. 

As a result we consider just that steps of splitting which deals with the remaining integral term. 
The corresponding equation reads 

f) r°° p -VR.\y\ 

-C(S,V R ,V L ,t) = -y/VRj o C(Sey,V R ,V L ,t)X R —^dy (22) 

for positive jumps and 

B , rO P -"L\y\ 



C(S,V R ,V L ,t) = -^/Vl [ C{Sey,V R ,V L ,t)\ L ^^-dy (23) 

J— oo \y\ 



for negative jumps. 

Making a change of variables x = log S and omitting dependence on dummy variables V R , Vl we 
can rewrite these two equations in a more standard form 

Q r°° e - v B,\v\ 

-C(x,t) = -V^^ C(x + y,t)\ R — TT ^dy (24) 

d_ , ., rrr- f° „, . ... e~ uM 

dt 



/-u e -"L\y\ 

C(x,t) = -VVl C(x + y,t)\ L — TT ^;dy 
J— oo \y\ 



To make it clear the above is not a system of equations but rather two different steps of the 
splitting procedure. 

Now an important note is that in accordance with the definition of these integrals we can rewrite 
the kernel as 

q r<x> e ~ v R\y\ 

-C(x,t) = -VVrJ o C(x + y,t)\ R 1+aR l y>0 dy (25) 

q r0 e -"L\y\ 

w C(S,t) = -VVl C{x + y,t)\ L —^-l y<a dy 
(ft j— oo \y\ 

This two equations are still PIDE or evolutionary integral equations. We want to apply our new 
method to transform them to a certain pseudo parabolic equations. 

First equation in the Eq. (1251) Assuming z = x + y we rewrite it in the form 

Q poo g — vr\z~ x\ 



-C(x,t) = -y/V R j C(z,t)X R 1+aR l z - x>0 dz (26) 
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In[3]:= 


Dl[p_] := 




^ p+1 




j Binomial [p + 1, i] k A (p + 1 - i) D [ x A p Exp [ -k x] HeavisideTheta [x] , {x, i}] 




Factorial [p] fz£ 


ln[5]:= 


FullSimplify [Dl [5] ] 


Out[5]= 


DiracDelta [x] 




FullSimplify [Dl [15] ] 


Out[6]= 


DiracDelta [x] 


Figure 1: Mathematica commands to verify the proposition |3.1 1 



At which Green's function is the kernel of the integral in the Eq. (|26p . i.e. 



To achieve our goal we have to solve the following problem. We need to find a differential operator 



(27) 



Al 



We prove the following proposition. 



e -v\y\ 



S(y) 



Proposition 3.1. Assume that in the Eq. a G I, and a < 0. Then the solution of the Eq. pJTU 

with respect to Ay is 



1 



Ap! 







9y 



l 



Xpl 



^ 1 dy l 



i=0 



■(1 + a) > 0, 



where Cf +1 are the binomial coefficients. 



Proof. As it will be shown later this result could be proven by taking Laplace transform of both 
parts of the Eq. (|27[) . It could be also verified using Mathematica commands given in Fig. Q] (we can 
check the above result for any positive integer p). □ 



Second equation in the Eq. ( 1251) For the second equation in the Eq. (j25j) it is possible to 
elaborate an analogous approach. Again assuming z = x + y we rewrite it in the form 



C(z,t)\ 



-u R \z-x\ 



\z — X 



l+a R 



-z-x<0 



dz 



(28) 



Now we need to find a differential operator A y which Green's function is the kernel of the integral 
in the Eq. (J2SD, i.e. 



a: 



We prove the following proposition. 



e 

\y\ 



(29) 



Proposition 3.2. Assume that in the Eq. i29\) a £ I, and a < 0. Then the solution of the Eq. h2§\ 

with respect to A~ is 



Xpl 







A, = -r-r \v - = — 



9y 



p+i 



Xpl 



■p+i 



0' 



^ ' 1 dy l 



i=0 



-(1 + a), 
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D2[p_] := 

Factorial [p] 

p+i 



V (-1) A i Binomial [p + 1, i] k A (p + 1 - i) D[ (-x) A pExp[-k (-x) ] HeavisideTheta [ -x] , {x, i} ] 



FullSimplify [D2 [5] ] 

DiracDelta [x] 



Figure 2: Mathematica commands to verify the proposition 13.21 



Proof. Using Laplace transform or Mathematica commands given in Fig. [2] we can check the above 
result for any positive integer p. □ 

To proceed we need to prove two other statements. 
Proposition 3.3. Let us denote the kernels as 

e -v R \z-x\ 

g + (z — x) = Xr- -j— l z - x>0 . (30) 

\z — x\ R 

Then 

A~g + {z — x) = 5(z — x). (31) 

Proof. 

a- +/ \ W d \ p+1 , . . i / d y +1 

A T g (z — x) = — — \ u — — g^(z — x) = —— v + — g^iz — x) 

x v ; Xp\ V dx J y ' Xp\\ d(z-x)J y ' 

= A+_ x g + {z-x) = 5{z-x) 
Proposition 3.4. Let us denote the kernels as 

-v h \z—x\ 



A z-x9 {z - x) = 5{z - x) 



□ 



g (z-x) = X L - t^-±z-x<q- (32) 

\z — x\ L 

Then 

Atg~{z-x) = 5{z-x). (33) 

Proof. 

,+ -, ^ 1 ( d \ P+l ^ 1 ( d Y +1 -i 
ATg (z — x) = — — \u + — g (z — x) = — — v — — g (z — x) 

x v ; Xp\ V dx) K ' Xp\\ d(z-x)J v ; 



□ 



A Itkin P Cafr s i n § pseudo-parabolic and fractional equations for option pricing in jump diffusion models 



10 



Transformation We now apply the operator A x to both parts of the Eq. ([25]) to obtain 

Q roc ( roc 

A~-^C{x,t) = -^/VrA~ J C(z,t)g + (z - x)dz = -^/VrIJ C(z,t)A~g + (z - x)dz + KJ 
= -y/^y^C(z,t)6(z-x)dz + n\ = ~VVRC(x,t) - ^JT R K 



(34) 



Here 



i=0 



z—x=0 



and a% are some constant coefficients. As from the definition in the Eq. (|30[) g(z — x) oc (z — x) p , the 
only term in the Eq. (|35p which does not vanish is that at i = p. Thus 

/ d p \ 

ll = V{x)[—g{z-x) \ z _ x=q = V(x)p\l {0) = 0; (36) 

With allowance for this expression from the Eq. (|34|) we obtain the following pseudo parabolic 
equation for C(x,t) 

A~^C(x,t) = -ly/vkC{x,t) (37) 

Applying the operator A% to both parts of the second equation in the Eq. (j28|) and doing in the 
same way as in the previous paragraph we obtain the following pseudo parabolic equation for C(x,t) 

At^-C{x,t) = - l -y/V L C(x,t) (38) 

4 Solution of the pseudo parabolic equation 

Assume that the inverse operator A~ l exists (see discussion later) we can represent, for instance, the 
Eq. (|37p in the form 

^C(x,t) = -BC(x,t), B=^^/Vr(A~)-\ (39) 
This equation can be formally solved analytically to give 

C(x,t) = e B{ - T ' t] C{x,T), (40) 

where T is the time to maturity and C(x,T) is payoff. Switching to a new variable r = T — t to go 
backward in time we rewrite the Eq. (|40|) as 

C{x,t) = e Br C(x,0), (41) 

Below we consider numerical methods which allow one to compute this operator exponent with a 
prescribed accuracy. First we consider a straightforward approach when a £ I. 
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4.1 Numerical method when a £ I 



Suppose that the whole time space is uniformly divided into N steps, so the time step 6 = T/N is 
known. Assuming that the solution at time step k, < k < N is known and we go backward in time, 
we could rewrite the Eq. (|40p in the form 

C k+1 (x) = e m C k (x), (42) 

where C k (x) = C(x, k6). To get representation of the rhs of the Eq. (|42|) with given order of approx- 
imation in 6, we can substitute the whole exponential operator with its Pade approximation of the 
corresponding order m. 

First, consider the case m = 1. A symmetric Pade approximation of the order (1,1) for the 
exponential operator is 

m 1 + B9/2 

1-B0/2 { ' 

Substituting this into the Eq. (|42|) and affecting both parts of the equation by the operator 1 — B8/2 
gives 

1 - ^Bo\ C k+l (x) = (l + ±B0\ C k (x). (44) 



This is a discrete equation which approximates the original solution given in the Eq. (|42p with the 
second order in 8. One can easily recognize in this scheme a famous Crank-Nicolson scheme. 

We do not want to invert the operator A% in order to compute the operator B because B is an 
integral operator. Therefore, we will apply the operator A% to the both sides of the Eq. (j44]). The 
resulting equation is a pure differential equation and reads 

[a--^9\ C k+ \x)=[A- + - 




Let us work with the operator A~ (for the operator A% all corresponding results can be obtained 
in a similar way). The operator A% contains derivatives in x up to the order p + 1. If one uses a 
finite difference representation of these derivatives the resulting matrix in the rhs of the Eq. (|45p is 
a band matrix. The number of diagonals in the matrix depends on the value of p = — (1 + or) > 0. 
For central difference approximation of derivatives of order d in x with the order of approximation q 
the matr ix will h ave a t least I = d + q diagonals, where it appears that d + q is necessarily an odd 
number dEberlvl (j2008T )b Therefore, if we consider a second order approximation in x, i.e. q = 2 in 



our case the number of diagonals is/=p + 3 = 2 — or. 

As the rhs matrix T> = A% — \/Vr9 / 4 is a band matrix the solution of the corresponding system of 
linear equations in the Eq. (|45p could be efficiently obtained using a modern technique (for instance, 
using a ScaLAPACK package) . The computational cost for the LU factorization of an N-by-N matrix 
with lower bandwidth P and upper bandwidth Q is 2NPQ (this is an upper bound) and storage-wise 
- N(P + Q). So in our case of the symmetric matrix the cost is (1 — an) 2 N/2 performance-wise and 
iV(l — aji) storage-wise. This means that the complexity of our algorithm is still O(N) while the 
constant (1 — or) 2 /2 could be large. 

A typical example could be if we solve our PDE using an x-grid with 300 nodes, so N = 300. 
Suppose apt = —10. Then the complexity of the algorithm is 60iV = 18000. Compare this with the 
FFT algorithm complexity which is (34/9)2AT log 2 (2A^) w 20900 0, one can see that our algorithm is 
of the same speed as the FFT. 



1 We use 2N instead of N because in order to avoid undesirable wrap-round errors a common technique is to embed a 
discretization Toeplitz matrix into a circulant matrix. This requires to double the initial vector of unknowns 
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The case m = 2 could be achieved either using symmetric (2,2) or diagonal (1,2) Pade approxi- 
mations of the operator exponent. The (1,2) Pade approximation reads 



1 + B9/3 



" 1 - 2B9/3 + B 2 9 2 /6 : 
and the corresponding finite difference scheme for the solution of the Eq. (|42p is 



(A-) 2 -lVvROA- + ±v R e' 



C k+1 {x)=A- 



C k {x). 



which is of the third order in 9. The (2,2) Pade approximation is 

Be = l + B9/2 + B 2 9 2 /l2 

l- Be/2 + B 2 e 2 /\2' 

and the corresponding finite difference scheme for the solution of the Eq. (|42|) is 

C k+1 (x) 



{A-? + \^/v R eA- + ^y R e< 



C k (x) 



(46) 



(47) 



(48) 



(49) 



which is of the fourth order in 6. 

Matrix of the operator (A x ) 2 has 21 — 1 diagonals, where I is the number of diagonals of the matrix 
A~ . Thus, the finite difference equations Eq. (|4"T|) and Eq. (|4^|) still have band matrices and could be 
efficiently solved using an appropriate technique. 



4.2 Stability analysis 

Stability analysis of the derived finite difference schemes could be provided using a standard von- 
Neumann method. Suppose that operator A~ has eigenvalues £ which belong to continuous spec- 
trum. Any finite difference approximation of the operator A% - FD(A~) - transforms this continuous 
spectrum into some discrete spectrum, so we denote the eigenvalues of the discrete operator FD(A~) 
as Ci, i = 1, N, where iV is the total size of the finite difference grid. 

Now let us consider, for example, the Crank-Nicolson scheme given in the Eq. (|45p . It is stable if 
in some norm II • II 




It is easy to see that this inequality obeys when all eigenvalues of the operator A% are negative. 
However, based on the definition of this operator given in the Proposition 13 . 2 (. it is clear that the 
central finite difference approximation of the first derivative does not give rise to a full negative 
spectrum of eigenvalues of the operator FD(A~). So below we define a different approximation. 



Case an < 0. Therefore, in this case we will use a one-sided forward approximation of the first 
derivative which is a part of the operator f v R — -J^j . Define h = (x max — x m i n )/N to be the grid 
step in the x-direction, N is the total number of steps, x m i n and x max are the left and right boundaries 
of the grid. Also define c k = C k (xi). To make our method to be of the second order in x we use the 
following numerical approximation 

dC k (x) -C k +2 + AC k +1 -?,C k 2 

~^x~ ~ 2h +U[h) (51j 



A Itkin P Caff^S pseudo-parabolic and fractional equations for option pricing in jump diffusion models 



13 



Matrix of this discrete difference operator has the following form 



M f = — 
1 2h 



( 


-3 


4 


-1 





...0 \ 







-3 


4 


-1 


...0 










-3 


4 


...0 


V 





0... 








"3 / 



(52) 



All eigenvalues of Mf are equal to —3/(2h). 

To get a power of the matrix M we use its spectral decomposition, i.e. we represent it in the 
form M = EDE', where D is a diagonal matrix of eigenvalues di,i = 1,N of the matrix M, and E 
is a matrix of eigenvectors of the matrix M. Then M p+1 = ED P+1 E', where the matrix D p+1 is a 
diagonal matrix with elements d P+1 ,i = l,N. Therefore, the eigenvalues of the matrix (vr — -§^) aR 
are [i/r + 3/(2h)] aR . And, consequently, the eigenvalues of the matrix B are 



Cb = y/VR\ R T(-a R ) {Wr + 3/(2h)] a * -!/««}. 
As an < and vr > it follows that Cb < 0. Rewriting the Eq. ([31]) in the form 



C k+ \x) 



i-\bo 



-1 



and taking into account that Cb < we arrive at the following result 

-i 



1 



1 



B9 



i + iw 



< 1. 



(53) 



(54) 



(55) 



We also obey the condition R (vr — ^) > 0. Thus, our numerical method is unconditionally 
stable. 



Case «l < 0. In this case we will use a one-sided backward approximation of the first derivative 
in the operator ( vl + -^J which reads 



dC k (x) 



+ 0(h I 



dx 2h 
Matrix of this discrete difference operator has the following form 



M b = — 
2h 



(56) 



/ 


3 








...0 \ 




-4 3 








...0 




1 -4 


3 





...0 


V 


0... 


1 


-4 


3 / 



(57) 



All eigenvalues of M b are equal to 3/(2h). Then doing in a similar way as above we can show that 
the eigenvalues of the operator B read 



Cb 



^LA L r(-a L ){K + 3/(2/ l )] c 



.OIL 



} 



(58) 



As ol < and vl > it follows that Cb < 0, R [vl + ^) > 0, and the numerical method in this 
case is unconditionally stable. 
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4.3 Numerical examples 



Here we describe two series of numerical experiments. In the first series we solve the equation 

Q r°° e ~ v R\v\ 

—C(x,t) = J^ C(x + y,T)\ R 1+aR dy, a R < -1 (59) 

by using FFT and the finite difference scheme constructed based on computation of the Eq. (|59p 
with qj;, £ I and interpolation as it was described in section |4~T1 

We solve an initial problem despite it is easy to consider a boundary problem as well. We consider 
a put option with time to maturity T = 30 days. As the terminal condition (we compute the solution 
backward in time) we chose a Black-Scholes put value at r = where the interest rate is r = 0.01, 
the volatility is 0.1 and the strike is K = 100. We create a uniform grid in time with Nt = 50 nodes, 
so 6 = T/Nt is the step in time. 



FFT. To apply an FFT approach we first select a domain in x space where the values of function 
C(x,t) are of our interest. Suppose this is x £ (— x*,x*). We define a uniform grid in this domain 
which contains N points: x\ = —x*,X2, ...xn-i,xn = such that X{ — X{-\ = h,i = 2...N. We then 
approximate the integral in the rhs of the Eq. (|59p with the first order of accuracy in h as 

roc p -vr\v\ h ^zi e - v nM 

/ C{x + y,T)\ R —^dy = -Y J Ci +j {r)f j , f^Xn—^+Oih 2 ). (60) 

Jo \y\ 2 jti-i \ x i\ 

This approximation means that we have to extend our computational domain to the left up to 
xi_tv = x\ — hN. 

The matrix |/| is a Toeplitz matrix. Using FFT directly to compute a matrix- vector product 
in the Eq. (|60p will produce a wrap-round error that significantly lowers the accuracy. Therefore a 
standard technique is to embed this Toeplitz matrix into a circulant matrix T which is defined as 
follows. The first row of F is 

Fl = (/0)/l)---)/jV-l,0,/i_jv, •••,/-!), 



and others are generated by permutation (see, for instance, Zhang and Wang ( 20091 )). We also define 
a vector 



C=[C 1 (t),...C n (t), 0^,0] 

N 



T 



Then the matrix-vector product in the rhs Eq. (160p is given by the first N rows in the vector V = 
ifft(fft(Fi) * fft(C)), where fft and ifft are the forward and inverse discrete Fourier transforms as 
they are defined, say in Matlab. In practice, an error at edge points close to x\ and xjy is higher, 
therefore it is useful first to add some points left to x\ and right to xjq and then apply the above 
described algorithm to compute the integral. We investigated some test problems, for instance, where 
the function C was chosen as C(x) = x so the integral can be computed analytically. Based on the 
obtained results we found that it is useful to extend the computational domain adding N/2 points 
left to x\ and right to x^ that provides an accurate solution in the domain x\,...,xn- The drawback 
of this is that the resulting circulant matrix has 4N x 4N elements that increases the computational 
work by 4 times (4ZV log 2 (4iV) 4(iV log 2 N)). 

In our calculations we used = 20, h = 2x*/iV regardless of the value of N which varies in the 
experiments. Then we extended the domain to x\ = —x* — h(N/2 — l),xjy = x* + h(N/2 + 1), and 
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Time to compute the solution 
FD N = 0.01 5547 sec, N = 256, alpha R = -1 .000000 
FD N /FFT N =0.077030 
FD N /FFT 2N =0.076615 
FD N /FFT 4N =0.073234 
FD N /FFT BN =0.064363 
FD N /FFT ]6N =0.044886 



I 

M/A 






FD-FFT^ 

FD - FFT 2N 

FD - FFT 4N 

FD-FFT BN 

FD-FFT 1BN 







-20 -15 -10 -5 5 10 



Figure 3: Difference (FD-FFT) in solutions of the Eq. ( 15T?|) as a function of x obtained using our finite- 
difference method (FD) and an explicit Euler scheme in time where the jump integral is computed using 
FFT. a R = -1. 

so this doubles the originally chosen value of N, i.e. N new = 2N. But the final results were analyzed 
at the domain x E (— »*, x*). 

Integrating the Eq. (|59|) in time we use an explicit Euler scheme of the first order which is pretty 
fast. This is done in order to provide the worst case scenario for the below FD scheme. Thus, if our 
FD scheme is comparable in speed with FFT in this situation it will even better if some other more 
accurate integration schemes are applied together with the FFT. 

FD. We build a fixed grid in the x space by choosing S m i n = 10 -8 , S max = 500, x\ = log(S m j n ), x^ = 
\og(S max ), h = (xn — x{)/N, N = 256. A one-sided forward approximation of the first derivative was 
used as it is defined in the Eq. (|51|) to approximate the operator in the Eq. (|59p . In the particular case 
considered here in our experiments Vr = 1, and the compensators in the Eq. (|21J) are not considered, 
because they could be integrated out at aR < and added to the diffusion terms. The Crank-Nicolson 
scheme Eq. (j54"j) was applied to integrate the Eq. ([59]) in time. 

Results The first series of tests was provided when an £ I and vr = 1,A_r = 0.2. The results of 
this series are presented in Fig. [3][6j 

In case otR = —1 in Fig. [3] the FFT solution computed with N = 256 provides a relatively big error 
which disappears with N increasing. It is clear, because the Crank-Nicolson scheme is of the second 
order in h while the approximation Eq. (|60j) of the integral is of the first order in h. Numerical values 
of the corresponding steps in the described experiments are given in Tab. [TJ 

Therefore, h? FD ~ 1ifft w - Actually, the difference between the FD solution with Nfd = 256 and 
the FFT one with N = ANpo is almost negligible. However, the FD solution is computed almost 13 
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Time to compute the solution 
FD N = 0.01 5501 sec, N = 256, alpha R = -2.000000 
FD N /FFT N =0.077568 
FD N /FFT 2N =0.077975 
FD N /FFT 4N =0.072487 
FD N /FFT BN =0.064698 
FD N /FFT ]6N =0.044509 



- 












I 






FD-FFT^ 

FD-FFT 2N 

FD-FFT 4N 

FD-FFT BN 

FD - FFT , 6N 















-20 -15 -10 -5 5 10 



Figure 4: Same as in Fig. El oir = —2. 



Time to compute the solution 
FD N = 0.023826 sec, N = 256, alpha R = -5.000000 
FDJFFT ==0.114911 
FD N /FFT 3N =0.119001 
FD^FT ==0.110342 
FD N /FFT BN =0.09661 1 
FD N /FFT 16N =0.069316 




-0.4 - 
-0.6 - 
-0.8 - 



-1.2 - 



-20 -15 -10 -5 5 10 



Figure 5: Same as in Fig. [3] an = — 5. 



times faster. Even the FFT solution with N = Nfd is 10 times slower than the FD on^l 



2 It actually uses 47V points as it was already discussed 
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Time to compute the solution 
FD N = 0.026058 sec, N = 256, alpha R = -6.000000 
FD N /FFT N =0.124215 
FD N /FFT 2N =0.127160 
FD N /FFT dN =0.121105 
FD N /FFT BN =0.109183 
FD N /FFT ]6N =0.075362 





FD-FFT N 

— FD-FFT 2N 
FD-FFT 4N 

FD - FFT aN 
FD-FFT 16N 









-20 -15 -10 -5 5 10 



Figure 6: Same as in Fig. El oir = —6. 





FD 2 56 


FFT 2 56 


FFT512 


FFT1024 


FFT 2 o48 


FFT 40 96 


h 


0.096 


0.1563 


0.078 


0.039 


0.0195 


0.00977 



Table 1: Grid steps h used in the numerical experiments 



For or = —2 in Fig. 2] we see almost the same picture. For or = —5 speed characteristics of 
both solutions are almost same while the accuracy of the FD solution decreases. This is especially 
pronounced for atR = — 6 in Fig. [6] at low values of x. The problem is that when cxr decreases the 
eigenvalues of matrix B in the Eq. (I44p grow significantly (in our tests at cxr = — 6 the eigenvalues 
are of order of 10 7 ), so in the Eq. (|55p the norm of matrix is very close to 1. Thus the FD method 
becomes just an A-stable. However, a significant difference is observed mostly at very low values of 
x which correspond to the spot price S = exp(x) close to zero. For a boundary problem this effect is 
partly dumped by the boundary condition at the low end of the domain. 

The second series of tests deals with otR G M using the same parameters vr = 1,Xr = 0.2. The 
results of this series are presented in Fig. [TlfTO Four point cubic interpolation is used to compute the 
value of C(x,t) at real or using the closest four integer values of or. 

It is seen that cubic interpolation provides pretty good approximation to the solution which is 
comparable with the FFT method in the accuracy and is faster in speed. Again, as we already 
discussed at or < 5 the accuracy of the FD scheme drops down even for ocr £ K, therefore the same 
picture is observed for atR € BL 

At — 1 < atR < (see Fig. [TT|) the difference between FD and FFT solutions surprisingly increases 
with N, used in the FFT method, increasing. To better understand what is the reason of that we 
fulfilled a test calculation of the integral in the rhs of the Eq. (|59p when C(x, r) is a known function, 
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Time to compute the solution 
FD N = 0.038135 sec, N = 256, alpha R = -1 .500000 
FD N /FFT N =0. 193368 
FD N /FFT 2N =0.190916 
FD N /FFT JN =0.179190 
FD N /FFT gN =0.1 59630 
FD N /FFT 16fJ =0. 104849 



_ FD-FFT N 

- FD-FFT 2N 

- FD-FFT^ 
_ FD-FFT BN 
- FD-FFT 




Figure 7: Difference (FD-FFT) in solutions of the Eq. f[5"9"j) as a function of x at or £ K obtained using 
our finite-difference method (FD) and interpolation and an explicit Euler scheme in time where the jump 
integral is computed using FFT. cur = —1.5. 



FD N /FFT N =0. 199443 
FD N /FFT 2N =0.200650 
FD N /FFT 4N =0.1 88271 
FD N /FFT gN =0.1 65543 
FD./FFT 1 „.=0.110114 




Figure 8: Same as in Fig. [71 an = —2.5. 
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Time to compute the solution 
FD N = 0.040235 sec, N = 256, alpha R = -3.500000 
FD N /FFT N =0.191160 
FD N /FFT 2N =0.192177 
FD N /FFT JN =0.1 82708 
FD N /FFT gN =0.161503 
FD N /FFT 1gfJ =0.1 05971 



_ FD-FFT N 
_ FD-FFT 2K 
_ FD-FFT 4K 
_ FD-FFT 8K 
_ FD-FFT,„ 




Figure 9: Same as in Fig. [7J or = —3.5. 



Time to compute the solution 
FD N = 0.041590 sec, N = 256, alpha R = -5.500000 
FD N /FFT N =0.201605 
FD N /FFT 2N =0.1 93793 
FD N /FFT 4N =0.192164 
FD N /FFT gN =0.163117 
FD./FFT 1R ..=0.1 10791 




FD-FFT N 
FD-FFT 2r 
FD-FFT 
FD-FFT ef 
FD-FFT 



Figure 10: Same as in Fig. [7J a R = —5.5. 



namely C(x, r) 



x. In this case this integral can be computed analytically which gives 

{x + y) 1+ dy = {xv R - a R )v^ R T(-a R ). 
o \y\ 



(61) 
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Time to compute the solution 
FD N = 0.038226 sec, N = 256, alpha fl = -0.500000 
FD N /FFT N =0.191429 
FD N /FFT 2N =0.188312 
FD N /FFT 4N =0.173474 
FD N /FFT gN =0. 152407 
F D /F FT = . 1 4490 




Figure 11: Same as in Fig. [3 or = —0.5. 



Then we apply the above described FFT approach and compare the numerical solution with the 
analytical one. The results of this test are given in Fig. [12j It is seen that FFT algorithm used in 
our calculations doesn't provide a good approximation to the analytical solutions at low N. So we 
expect this behavior of the FFT method occurred in our numerical experiments at ocr = —0.5, but 
this doesn't explain the observed effect. 

A plausible explanat i on is t hat at cxr close to the integral kernel becomes singular. That is why 



Cont and Voltchkoval (120031 ) the part of the infinitesimal generator corresponding to small jumps 



m 

is approximated by a differential operator of second order (additional diffusion component). As we 
didn't use this technique here, an increase of N forces the distance between y = and the closest 
FFT node boundary to become smaller, thus the kernel becomes larger. 

The other reason for the FD solution to differ from the FFT solution is that at — 1 < cir < we 
don't use the option values computed at an = (remember, this is a special case that was discussed 
earlier). Thus, instead of interpolation we use extrapolation that certainly decreases the accuracy of 
the FD solution. We will resolve this problem in the next section. 

At the end of this section we present the option values computed using such a scheme as a function 
of x obtained in the same test (Fig. [T3|) . 



5 General case 

If we take a more close look at the propositions 13.11 and 13.21 we could recognize that the assumption 
a £ I could be neglected while both propositions will remain valid. This could be easily seen based 
on the following equalities 



Proposition 5.1. Assume that in the Eq. |^7[ ) a 6 R, a < —1. Then the solution of the Eq. \27^ 
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FFT computation of a test integral 
Depicted: Numerical solution - exact solution 
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Figure 12: FFT computation of a test integral in the Eq. (I6T|) 



Option value at various ci R and lime to compute it 



FD 


= 0.039698 SE 




-4.500000 


FD 


= 0.006165 SE 




-4.000000 


FD 


= 0.035964 SE 


c, alpha R = 


-3.500000 


FD 


= 0.005952 SE 


c, alpha R = 


-3.000000 


FD 


= 0.036014 SE 


c, alpha R = 


-2.500000 


FD 


= 0.005499 SE 


c, alpha R = 


-2.000000 


FD 


= 0.034729 SE 




-1 .500000 


FD 


= 0.005699 SE 


c, alpha = 


-1 .000000 


FD 


= 0.033833 SE 




-0.500000 



Figure 13: Option values computed using such a scheme as a function of x obtained in the same test 



with respect to A£ is 
1 



XT( P + 1) 



d_ 

dx 



p+i 



1 



XT(p + l) 



dx 1 



P 



-(1 + a) > 0, 
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where Cf +1 are the generalized binomial coefficients which could be exp ressed via Gamma fun c tion, and 
fractional derivatives are understood in the Riemann-Liouville sense /(Oldham and Spanien If 1974 )) 

Proof 1. Taking Laplace Transform of the expression A + f(x) we obtain 



0' 



ErjP+ 1 u P+ 1 - i __ 
dx l 



i=0 



i 



XT(p + 1) 



p+l-i J 



i=0 



1 



XT(p+l 



{v + sT +1 £ s f(3 



Now, as 



/(*) = A: 



-u\x\ 



il+^- L 3/>0 



and 

we obtain 

And thus A+f(x) = 5(x) 



cA x. 



-v\x\ 



|l+a 



lx>o^ = XT(p+l)(v + s) 



-(i+p) 



C s {At f{x)) = 1 = C s 5(x) 



For the operator A x the proof is similar. 
Another proof is based on a different idea. 

Proof 2. As it is well known a shift operator in L2 space could be represented as follows 

6. ' 8 



SO 



exp a— 

ox 



6 a f(x) = f(x + a). 
Therefore, the integrals in the Eq. (|24p could be formally rewritten as 

AiC(x,t), Ai = 
A 2 C(x,t), A 2 = 



/•oo 

/ Xr 


e 


-"ii\y\ 


exp 


( d 

\ V dx 


y 




1 — oo 


e~ 


-vl\v\ 


exp 


( y— 

^ dx 


\y 





We can compute these integrals assuming that d/dx is a constant. This gives 



Ai = X R r(-a R ) v R 



A 2 = X L T(-a L ) ( v L + — 



dx y 



(a) < 0,R(v R -d/dx) > 
a) < 0,M(z^ L + > 0, 



Csf(x) 



□ 



(62) 



(63) 



(64) 



where under a real part of differential operator we will understand the real part of the maximum 
eigenvalue of finite difference matrix which approximates this differential operator (see below). 
A simple observation shows that 



A! = {A-)-\ A 2 = {A+) \ 



which finalizes the proof. 



□ 



A Itkin P Caff^S pseudo-parabolic and fractional equations for option pricing in jump diffusion models 



23 



This means that the whole analysis of the previous sections made in the case a G I is still valid 
for arbitrary a £ K, a < 0. Moreover, we could now extend this proof for the whole range of a < 2. 
In order to do that we have to consider the whole integrals in the Eq. (I2ip . This is because in the 
case of jumps with infinite activity or infinite variation the second and third integrands can not be 
integrated out, because they do not exist. 

If we apply the second transformation to the first equation in the Eq. (|24|) the result is given by 
the following proposition. 

Proposition 5.2. The PIDE 



^-C(x,V r ,V l ,t) 







VVr 

is equivalent to PDE 

VrXrT(-(xr 







C(x + y, Vr, V l ,t) - C(x, Vr, V l , t) - —C(x, Vr, V L) t)((* - 1) 



e -"R\v\ 



\y 



—C(x,Vr,V l ,t) 



VR 



M(a R ) < 2, R(u R - d/dx) > 0, R(v R ) > 1 
In special cases this equation changes to 



q- x ) - »r + k r - - i) Qfl ] q- x ) c{x > Vr ' Vl > t) > 

(66) 



d_ 

dr 



and 



C(x,V r ,V l ,t) 



VrXr < log(v R ) - log v R 



Ot 



a R = 0,R{vR-d/dx) > 0, 



C(x,V R ,V L ,r) 



d_ 

dx 

> 1, 



+ log 



VR 



d_ 

dx 



C(x,V R ,V L ,r) (67) 



VrXr\ ~ Vr log Vr + (vr 





d_ 

dx 



log Vr 



d_ 

dx 



(68) 



+ [vRlogVR - (V R - l)l0g(^ - 1)] —jC(x,V R ,V L ,T) 

a R = l,R(d/dx) < 0,R(vr) > 1, 
where logarithm of the differential operator is defined in a sense of (Bakas et al. (199&) ). 
Proof. We again use the shift operator introduced in the Eq. (|62p to rewrite the Eq. (|65p as 
d 



dr 



C(x, V R , V l , t) = BtC(x, V R , V L , r) 



Bi = VVr 



-vn\y\ 



-dy 



(69) 



Formal integration could be fulfilled if we treat a differential operator ^ as a parameter. As it 



could be verified the result is that given in the Eq. (|66j) . Same method is used to prove the formulae 
given in the special cases ocr = and ocr = 1. □ 



Also notice that at or = from the very beginning the last term in the Eq. (J65J) can be moved 
from the integral to the diffusion part of the Eq. (|2ip because the remaining kernel converges at y = 0. 
If we do so, at this special case the integrated equation transforms to 



d_ 

dr 



C(x,Vr,V l ,t) 



VrXr < log(v R ) - log v R 



d_ 

dx 



C(x,Vr,V l ,t) 



(70) 



a R = 0,R(v R - d/dx) > 0,R(i/ fl ) > 0, 
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This form is more useful as we show later when elaborating a numerical method to solve it. 



The same approach could be utilized for the second equation in the Eq. (|24p . and the result is 
given by the following proposition. 

Proposition 5.3. The PIDE 





—C(x,Vr,V l ,t) 







C{x + y, V R , V l ,t)- C(x, Vr, V l ,t) - —C(x, V R , V L ,r)(e y - 1) 



Al 



-vl\v\ 



1 2/1 



dy (71) 



is equivalent to PDE 



^-C(x, V R , V l ,t) = VV[X L T(-a L ) Mv L + JL\ * - + [v a L L - (»l + l)° L } ^ } C(x, V R , V L , r), 



R(a L ) < 2, R{u L + d/dx) > 0, R{u L ) > 0. 
In special cases this equation changes to 

^-C(x,Vr,V l ,t) = -y/V L \ L | log (v L + J;) " log(^) " log 
a L = 0, R{u L + d/dx) > 0, R{u L ) > 0, 



vl + 1 



d_ 

dx 



(72) 



(73) 



and 



d_ 
dr 



C(x,Vr,V l ,t) = yJV~ L \ L {- v L \ogv L 

d d ( d \ 

+ K V L ~ {VL + 1) log(^L + 1 J\fo + fa + ^ l0g Y L+ dx~) l C ^' VL,T ^ 

•x) < 0, R(u L ) > 0, 



(74) 



where logarithm of the differential operator is defined in a sense of (Bakas et al. (199&) ). 
Proof. The proof is similar to that given in the Proposition 15.21 



□ 



Again at ar, = we can move out the last term in the Eq. (|65p from the integral to the diffusion 
part of the Eq. (I21D because the remaining kernel converges at y = 0. If we do so, at this special case 
the integrated equation transforms to 



d_ 

dr 



C(x,Vr,V l ,t) 



V l \ l I log (v L + ^\ - \og{u L ) 
a L = 0, R(u L + d/dx) > 0, R(u L ) > 0, 



(75) 



We will use this form later when elaborating a numerical method to solve this equation. 

It is important to underline that the integration in the Proposition 15.21 for positive jumps could be 
done if M(^i?) > 1 while in the Proposition 15.31 for negative jumps - if R(vl) > 0. In the special cases 
ctR = 1 this limit could be extended to R(z<r) > 0, however it gives rise to a complex values of the 
coefficients in the rhs of the Eq. (|74|) . Therefore, we keep the above constraint M(z^r) > 1 unchanged 
in this case as well. 
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Similar representations were obtai ned first in iBovarchenko and Levendorskii (|2002h and later in 
Cartea and del Castillo-Negrete ( 20071 ) using a characteristic function approach. For instance, the lat- 



ter authors considered several Levy processes with known characteristic function, namely LS, CGMY 
or KoBoL. Then using Fourier transform they managed to convert the governing PIDE (same type 
as the Eq. (|21j) but for the Black-Scholes model with jumps) to a fractional PDE. In their notation 
our operator Ai is represented as 



Ai oc (- 



S DS* ( e -"*C(x,t)), 



and operator A2 as 



A 



2 oc e VL ooK L 



-''L 



C(x,t)) 



(76) 



(77) 



So to compare we have to note that aside of the different method of how to derive these equations 
our main contribution in this paper is: 



1. 



5. 



Sp ecial cases a r = 0, 1,cy/ = 0, 1 a re no t considered in ICartea and del Castillo-Negretd (|2007l ) . 
In Boyarchenko and Levendorskii ( 20021 ) a corresponding characteristic function of the K0B0L 
process was obtained in all cases for a < 1. However, the authors did not consider numerical 
solution of the fractional PDE. In this paper we derive a fractional PDE for all a < 2 and 
propose a numerical method for their solution. 

We proposed the idea of solving FPDE with real an < 0, oll < by using interpolation between 
option prices computed for the closest integer values of an,ai. For the latter we first used to 
transform the fractional equation into a pseudo-parabolic equation. Then for the solution of this 
PPDE an efficient FD scheme is constructed that results in LU factorization of the band matrix. 

Also jumps up and down are c onsidered separately so the mode l in us e (SSM) is slightly different 
from the model considered in lCartea and del Castillo-Negretd (120071 ). 



In lCartea and del Castillo-Negretd (120071 ) a Crank-Nicolson type numerical scheme was proposed 
to solve the obtained FPDE in time while discretization in space was done using the Grunwald- 
Letnikov approximation which is of the first order in space. Here for fractional equations with 
2 > an > 0, 2 > «i > we obtain the solution using our new scheme which preserves the second 
order approximation in time and space. 



As it is known from rece n t papers (|Abu-Saman and Assail (120071 ). iMeerschaert and Tadjeran 
(j2004l . \20Q(h . ISousal (120081 ). iTadjeran et al.1 (|2006l )). a standard Grunwald-Letnikov approxima- 
tion leads to unconditionally unstable schemes. To improve this a shifted Grunwald-Letnikov 
approximation was proposed which allows construction of the unconditionally stable scheme of 
the first order in space. H Here we use a different approach to derive the unconditionally stable 
scheme of higher order. 

6. We show that when considering jumps with finite activity and finite variation despite it is a 
common practice to integrate out all Levy compensators in the Eq. (|21[) in the integral terms 
this breaks the stability of the scheme at least for the fractional PDE. Therefore, in order to 
construct the unconditionally stable scheme one must keep some other terms under the integrals. 
To resolve this in Cartea (2007) the authors were compelled to change their definition of the 
fractional derivative (see below). 

7. Our approach could be easily generalized for a time-dependent Levy density. 



3 A second order approximation could in principle be constructed as well, however resulting in a massive calculation of 
the coefficients. That probably stopped the scientists to further elaborate this approach. 
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6 Numerical method 



Let us consider a general case which is given by the Eq. (|66p and Eq. (|72p □. We first discuss how to 
construct an unconditionally stable scheme of the second order in space and second or higher order in 
time. Then we consider some peculiarities of implementation of the derived finite difference schemes. 

6.1 Case ocr = or oll = 0. 

This extreme case corresponds to the familiar Variance Gamma model. In this case the integrals 
in the Eq. ([65p and Eq. (|71|) exist if we keep just first two terms under the integral. Therefore we 
could integrate out the last term R oc J|C(x,r)(e y — 1). This term then will become a part of the 
convection part of the total PIDE and therefore we will not consider it here, assuming that we use a 
splitting technique and know how to solve the remaining convection-diffusion equation. 
Then the Eq. (|66p could be written in the form of the Eq. (j39|) with 



In practical computation of the rhs operators we exploit a modification of our interpolation method 
which was described above. First, note that typical values of \r, \l as well as Vr, Vl are limited, i.e. 
normally \r < M, Al < M, Vr < M, Vl < M where M could be chosen in the range, say 3-5. Second, 
if we solve a general jump-diffusion equation using some kind of splitting methods, the time step of 
integration 9 in the Eq. (|79p is determined by the time step used at the integration of the diffusion 
part. This means that 6 is usually small. Therefore, it is pretty reasonable to assume that in the 
Eq. (|79p m < 2. Next, as follows from the definition of the fractional derivatives, the operators in 
the Eq. (|79p are continuous in m. Therefore, we could solve the Eq. (|79p for m = 0, 1, 2 and then use 
quadratic interpolation to get the solution given the real value of m, and the condition m < 2. Note, 
that m = is a trivial case so the solution C k+1 {X) = C k (x) is already known. 

Note a choice of m = —1. On the one hand this is very attractive because then the solution of 
the Eq. (|79p is already found. On the other hand at m < the scheme in the Eq. (|79p becomes 
explicit which breaks its unconditional stability. Apparently the best one can achieve in this case is 
to use a central difference approximation for the first derivative. Then it is possible to show that all 
eigenvalues of the rhs matrix have their real value equal to one. Thus the stability of the scheme is 
questionable. 

We now construct a stable FD scheme to solve the first equation in the Eq. (I79p . Similar to 
what was already discussed in the previous section a forward second order approximation of the first 

4 In principal one can eliminate special cases when one of the following conditions is valid o.r = 0, = 0, ccr = 1,ul = 1, 
by just substituting, say an = e << 1 instead of olr ~ 0, an = 1 + e instead of an = 1 etc 




(78) 



Therefore, integrating it we obtain an explicit form of the Eq. (I42p 




(79) 



A Itkin P Caff^S pseudo-parabolic and fractional equations for option pricing in jump diffusion models 



27 



( 1 A \ _m 

derivative has to be chosen. Then the eigenvalues of the discrete operator I 1 — jt^j^:) are 



1 + — . (80) 

2hu R J v ; 

We need to guarantee that || (1 — jt^-^J II < 1- Thus, if vr < 1 this FD scheme is stable 
at h < 3/ [2(1 — vr)], and if v R > 1 - it is unconditionally stable. As follows from the Proposition 
Eq. (|5.2p M(i/r) > 1, therefore the scheme is unconditionally stable. 

After this discretization the matrix of the lhs operator becomes one-sided tridiagonal if m = 1, 
and one-sided pentadiagonal if m = 2. Therefore this equation can be efficiently solved with the total 
complexity 0(N(2m + 1)). 

To preserve monotonicity of the solution for the second equation in the Eq. ([79]) a backward second 
order approximation of the first derivative has to be chosen. This approximation was also already 
introduced in the previous section. Then C k+l (x, m) can be computed as a product A m ■ C k (x), where 
A m is a band matrix with 2m + 1 diagonals. So the complexity of this is also 0(N(2m + 1). 

Based on these results we extend our numerical test described in the previous section to the case 
oir • = 0. However, to preserve convergence of the integral now instead of the Eq. f)59[) we have to use 
the extended equation 



q r°° e -VR\y\ 

—C(x,t) = J [C(x + y,r)-C(x,T)}\R 1+ctR dy, a R < -1 (81) 



q r°° e -"ii\y\ 

\y\ 

We again compare the FFT solution of the Eq. (|8ip with that obtained based on our method. 



FFT. It should be underlined that the presented simple FFT algorithm completely loses its accuracy 
when or — > 0. Therefore, instead of or = we will chose real or = —0.5. We again define a uniform 
grid in the domain (— x*,x*) which contains points: x\ = — x*,X2, ...xn-i,xn = x* such that 
Xi — Xi-i = h,i = 2...N. We then approximate the integral in the rhs of the Eq. (|8ip with the first 
order of accuracy in h as 

poo e -va\v\ N ~ % 

/ [C(x + y, t) - C(x, t)] Xr 1+oir dy = h ^ C i+j (T)fj - C(x, T)\ R v% R T(-aR), 
Jo \y\ j=i-i 



-VR\xj\ 



U - ^-nr^ + °(^) ( 82 ) 

\Xj\ 

The matrix-vector product in the lhs of the Eq. (|82p is computed using FFT as it was described 
in the previous section. 

FD. We solve the Eq. (|8T|) using interpolation in or between the points or = 0,-1, —2, —3. At 
or = we use the FD scheme in the Eq. (|79p . At aR < we again use our approach of construction 
of the pseudo-parabolic equations (see propositions 13.31 13.4p . and instead of the Eq. (I39p now obtain 

^-C(x,t)=BC(x,t), B=^{A-)~ l -X R u^T(-a R ). (83) 



Further we use the Crank-Nicolson scheme Eq. (|44p which now reads 

A- \o\ C k +\x) = (l- ±\ RV %T(-a R )0 A- + \o\ C k {x). (84) 



1 + 7j A fl i#T(-afl) 
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Time to compute the solution 
FD N = 0.037912 sec, N = 256, alpha fl = -0.500000 
FD N /FFT N =0.179801 
FD^/FFT 2N =0.185900 
FD N /FFT 4N =0.179917 
FD^/FFT aN =0.152166 
FD N /FFT ]6N =0.104247 




Figure 14: Difference (FD-FFT) in solutions of the Eq. (I6TJ|) obtained using our finite-difference method 
(FD) and an explicit Euler scheme in time where the jump integral is computed using FFT. ctR = —0.5. 



The stability analysis could be provided similar to what we did in the previous sections. Again it 
is easy to show that the forward one-sided approximation of the operator A% given in the Eq. f)51 1) 
guarantees the unconditional stability of the above scheme. 



Comparison. The results of this test are given in Fig. Q3J This could be compared with the results 
presented in Fig. QTJ The difference is that now instead of extrapolation we use interpolation, because 
we are able to solve our test problem numerically at ctR = 0. Surprisingly the difference in the FFT 
and FD solutions slightly increases in case of interpolation. The FD solution is still faster than the 
FFT, and as follows from the above analysis - more accurate. 



6.2 Case oir = = 1. 

This is a case of jumps with infinite variation and infinite activity. Therefore we have to keep the 
whole integrals in the Eq. (|65p and Eq. (|7ip . i.e. in each integral we can not integrate the last term 
out because otherwise the integral does not converge. 

Let us remind that as follows from the Proposition 15.21 in this case the original PIDE Eq. (|65p is 
equivalent to the PIDE 



d_ 

dr 



C(x, V R , V l ,t) = t/VrXr^ - v R log vr + [y R - — ) log ( vr - — J 

d -i 

+ WRlogUR - (ur - l)log(z^ - 1)] —jC(x,V R ,V L ,T) 
R{d/dx) < 0,R(vr) > 1, 



(85) 
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while from Proposition 15.31 the PIDE Eq. (|74p is equivalent to the PIDE 

^C{x,Vr,V l ,t) = ^/V~ L \ L [-v L \ogv L (86) 

R(8/dx) < 0, R(y L ) > 0. 
For the following we need to prove the following Proposition. 
Proposition 6.1. The following identity holds 

-v R log v R + (y R — — ) log ( v R - — J + [i/r log i/fl - (v R - 1) log(zv R - 1)] — 

Proof. To prove this we one have to note that 

oo e -vn\y\ e -VR\y\ 



-dv 



and then use Proposition 15.21 with a R = 0. □ 

In a similar way we can prove the following proposition 
Proposition 6.2. 

d d ( 8 

-v L log v L + \v L log v L - (v L + 1) \og{v L + + K V L + -Tfo) log I v L + — 



y jlog(^i) - log + + log 



vl ) dx 



>dv (89) 



□ 

These two identities gives us an idea of how to construct a FD numerical method for solving the 
Eq. (|85[) and Eq. (|86p . First we rewrite the Eq. (|85p and Eq. (|86[) in the form 



^-C(x,V r ,V l ,t) =1L r C(x,V r ,V l ,t) (90) 
Lr = \/^?A jR y i log{u R ) - log f z/ B - \ + ^log ^— ^- \dv 



Ll = \/Vl^l y jlogK) - log ^ L + ^) + b g ( 



i/£, I dx 



We already know how to solve these equations if the operators Lr and Lx, do not contain the 
integrals. We want to utilize this approach by proceeding with the following steps. 
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Step 1. First we truncate the upper limit in the integral to some z/*. This could be done because 
the integral in the Eq. (|90p is well-defined and at vr — > oo the integral kernel tends to zero as 



(91) 



lim L r C(x,V r ,V l ,t) = VVrXr-^ (~ + ^ ) +0(1/4) 
vr^oo 2v R V ox ox 2 1 



At the interval (is, z/*) we approximate the integral in v using some quadrature formula, for instance, 
the well-known Simpson formula (higher-order approximations of even adaptive quadratures could 
definitely be used as well). So we partition the interval {y, i/*) into an even number of intervals M all 
of the same width h = (z/* — v)/M. Then operators in the Eq. (|90j) transform to 



M 



M 



(92) 



8=0 



i=0 



U,R = a,VV~nX R U -^f^o g (, h R) - log (v hR + (log ^-i) |- J 

Hl - «^A L ^|log(,, L ) - log (u itL + A) + log (^±1) |- j, 
Oi = l, i = 0,M, Oi = 2, i = 2,4...M-2, a* = 4, i = l,3...M-l. 

Step 2. Each operator in the Eq. (|92p is a sum of M operators which commute with each other. 
Therefore, the solution of the Eq. ()90j) reads 



C(x,Vr,V l ,t) = exp 
C(x,Vr,V l ,t) = exp 



■ Af 
.i=0 

.i=0 



C(x, V R , V l , 0) = J] e L ^ T C(x, Vk, Vl, 0) 

i=l 
A/ 

C(z, V R , V l , 0) = J] e L ^^C7(x, Vr, Vl, 0) 



(93) 



i=l 



Using a splitting technique (see, for instance, Lanser and Verwer ( 19981 ). Yoshida ( 199dl )) we can 
represent this equation in the form 



Ci(x, V R ,V L ,d) = e h ^ T C(x, V R , Vl, 0) 
C 2 (x, V R , V l , 9) = e^ T d{x, Vr, V l , 9) 

C M {x, Vr, V l , 9) = e h ^ T C M -i(x, Vr, V l ,0) 
C(x, V R , V l , 9) = C M (x, Vr, V l , 9) 



(94) 



and similarly for the operator Lj,. 



Step 3. Each equation in the Eq. (|94p is very similar to that corresponding to the case a = (see 
the previous section). The only difference is that the operators L^r now contain an extra term L3 i r = 

(log V%, y. ~r~ j ~Bxi anc ^ ^he operators Lj^ now contain an extra term IL^x, = (log ~§^- We 

can apply splitting to these operators similar to as we did in the above. Further by analogy with what 
was already discussed in the previous sections devoted to Pade approximations, these terms e h3 ' iM 9 
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and e L3 ' 4 - L # could be approximated with the second order of accuracy in 6 by using the Eq. (|43|) . 
Finally, each equation in the Eq. (|94p reads 



C k _f{x) = C k {x) 



(95) 



C?_i(x) 



1 



m; 



C k+1 (x) 



1 d 



vi,r dx 



Ci*(x), 



0,...,M, 



m = ai^V R \ R 



v* - v 
3M 



C\f\x) 



We can chose the number M to guarantee that the value of rrii is less than 2 and then use 
interpolation solving the above equations at rrii = 0, 1, 2. 

Similar scheme could be constructed for the operator which reads 



C^x) = C7 fc (x) 



(96) 



1 + fW? fe+1 

1 ^~^3,i, L# 



Cf +1 (x) 
C k+1 (x) 



1 + 



1 5 



fi,L dx 



C k +\x), i = 0,...,M, 



mi = ai\/V L \ L 



is* — v 



3M 



C\t\x) 



Step 4. To construct an unconditionally stable scheme in x we have to chose approximation for 
the first derivative in the Eq. (|95p . If we rewrite this equation in the form 



C^ix) = C k {x) 



(97) 



log 



Vi,R 



dx 



C 



k+l. 



1 + 



log 



Vi,R - 1 \ d 



Vi,R 



dx 



Cf_i(x), i = 0,...,M 



1 ^ 



C*(s) = C£(s) 



C k+1 (x) = CH\x) 



fe+i. 



it becomes obvious that the derivative in the second equation in the Eq. (|97p should be approximated 
by using a backward one-sided second order divided difference. For the derivative in the third equation 
one has to use a forward approximation. 

Similarly we rewrite the Eq. (|96|) in the form 



C k _f{x) = C k {x) 



(9£ 



log 



Vi.L + 1 \ 5 



1 + 



^i,L / dx 
1 5 



C£ +1 (x) 



1 + 



rrii 



log 



^i,L / dx 



Ct^x), i = 0,...,M 



vll dx 



Cr i (x) = CUx) 
C k+1 {x) = C k f\x) 



and use a forward approximation for the derivative in the second equation in the Eq. 
backward approximation in the third equation. 



and the 
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Solution as a function of v. 
Computational time = 0.613768 sec, N == 256, a = -0.500000, M = 30.000000 




'. = 300 
■. = 700 



Figure 15: Difference in solutions of the Eq. (I97p obtained at various and that at 
a R = 1. 



5 at M = 30 and 



The matrix in the rhs of the second equation in the Eq. (|98p is upper tridiagonal. The matrix in 
the rhs of the third equation in the Eq. (|98p is lower tridiagonal at mj = 1 and lower pentadiagonal 
at mi = 2. The total complexity of the algorithm as compared with the case a = is: one extra 
equation at each step, M steps instead of just one in the case a = 0. Therefore, using the results 
given in Fig. [J3] we can expect that at M = 30 this algorithm is about 3 times slower than the FFT. 
On the other hand it provides the second order approximation in both space and time, and does not 
require to re-interpolate the FFT results to the FD grid which was previously used to find solution 
for the diffusion part of the original PIDE. 

To verify this we provided two numerical experiments. In the first experiment varied while 
h = {y* — vr)/M was chosen to be constant. At = 5 we chose M = 30. The other parameters are 
same as in the previous numerical experiments reported in the above. This results are presented in 
Fig.dU 

The computational time rawly increases by the factor M/2, i.e. for M = 30 it is almost same as 
for the corresponding FFT. It is seen that an appropriate value of should be more than 300. 

In the second experiment we fixed the value = 300 and varied M to see at which M one could 
expect to get convergency. These results are presented in Fig [161 As it is seen M = 80 seems to be 
sufficient to obtain the convergency. The computational time in the case M = 81 is 1.4 sec which 
if compared with that given in the Fig Q3] is 3.6 times more than that for the FFT. Thus, in this 
case our algorithm is almost 4 times slower than the FFT. As it was already mentioned this could be 
compensated a) by the second order of accuracy in space and time, and b) no need for re-interpolation 
of the FFT results to the FD grid. On e more advantage is that we don't need to treat the point y = 
in a special way as it was done, say in ICont and Voltchkoval <|2003h . 

Note, that as we use M steps in the splitting scheme, the error in time becomes 0(M9 2 ) that 
could kill the second order of approximation. Therefore, for instance, in the Eq. (|95|) it is better to 
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Solution as a function of M 
Computational time = 0.368299 sec, N = 256, a R = 1 .000000, v, = 300.000000 




M = 21 

M - 81 

M = 101 

M = 201 

M - 401 



-20 -15 -10 -5 5 10 



Figure 16: Difference in solutions of the Eq. f!97|) obtained at various M and that at M = 21 at = 300 
and an — 1. 



use a third order approximation in time (see the Eq. (|46p ). Accordingly the second equation in the 
Eq. (|97l) will become 
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To preserve the third order of approximation in time the third equation in the Eq. (|97|) should 
now be solved at m = 0, 1, 2, 3 and then cubic interpolation to the actual value m 8 will give the final 
solution. This scheme increases the total computational time by about 10%, however the accuracy in 
time increases to 0(M6 3 ). 



7 Conclusion 

From the numerical point of view the proposed approach has an advantage as compared with the 
methods mentioned in the Introduction. Indeed, first we managed to reduce the original evolutionary 
integral equation to a pure differential equation. Second, this equation could be formally solved 
analytically. To compute the operator exponent we applied a Pade approximation technique. This 
eventually allowed us to derive finite difference equations which approximate the original solution 
with the necessary order. This equations could be solved at the same grid as the diffusion part of 
the original PIDE thus eliminating problems inherent to the FFT methods. In addition, despite 
the original integral term is non-local, the rhs matrix T> of the system of linear equations obtained 
by applying our approach is a band matrix in case of integer otR^, i.e. it corresponds to a local 
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approximation of the option price. Also we demonstrated that at a < the complexity of our 
algorithm is much lower than that of the FFT while the accuracy is much better. 

The complexity of the solution at a = 1 is higher than that of the FFT. This in part is compensated 
by few factors: our algorithm provides the second order approximation in both space and time, and 
it does not require to re-interpolate the FFT results to the FD grid which was previously used to find 
solution for the diffusion part of the original PIDE. 

Using this technique the solution at 2 > a > 1 could be obtained by using extrapolation given the 
solution at a = 1, 0, —1. 

It is interesting to know what are real values of a. In lBul (120071 ) the author used to calibrate the 
CGMY model to S&P 500 historical call option prices. The market prices were chosen from June 2007 
to December 2008. The strike is from 1300 to 2000 with the increment of 25 from 1300 to 1700 and the 
increment 100 from 1700 to 2000. The index closed price is 1536.34. The fou nd CGMY param eters 
were CGMY C = 0.0156, G = 0.0767, M = 7.5500, Y = 1.2996, i.e. a = 1.3. In lCarr et alT(|2005l l the 
option prices of S&P 500 were also calibrated using CGMY model which gave the values of a in the 
range (-0.39,-0.42). 
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